Z 2 -vortex lattice in the ground state of the triangular Kitaev-Heisenberg model 
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Investigating the classical Kitaev-Heisenberg Hamiltonian on a triangular lattice, we establish 
the presence of an incommensurate non-coplanar magnetic phase, which is identified as a lattice 
of Z2 vortices. The vortices, topological point defects in the SO(3) order parameter of the nearby 
Heisenberg antiferromagnet, are not thermally excited but due to the spin-orbit coupling and arise 
at temperature T — > 0. This Z2-vortex lattice is stable in a parameter regime relevant to iridates. 
We show that in the other, strongly anisotropic, limit a robust nematic phase emerges. 
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Mott insulators with orbital degrees of freedom [1] have 
long been a topic in condensed-matter research, because 
the physics of various transition-metal oxides is governed 
by anisotropics inherent in electronic orbitals. Apart 
from the drive to understand these materials, the lowered 
symmetry can lead to new physical phases, e.g., (spin-) 
orbital liquids [2-5], topological Mott insulators [6], or 
fractional Chern insulators [7] . Orbital anisotropics have 
additionally attracted interest as a highly desirable fea- 
ture in potential realizations of protected qbits [8, 9]. 

Iridates have been pointed out as a possible realiza- 
tion of strongly direction-dependent orbital Hamiltoni- 
ans [10], because their 5d electrons feel both substan- 
tial correlations and strong spin-orbit coupling. The 
spin-orbit coupling makes the total angular momentum 
the relevant quantum number to describe the low-energy 
states with a single hole in the 5d-t2 g levels [10-12]. The 
j = 1/2 of the low-energy states takes the role of a spin 
(and will be referred to as such from now on), but the 
orbital component heavily influences couplings between 
'spins' [10, 13]. While interactions in the Mott-insulating 
square planes of Sr2lr04 are isotropic, they are highly 
anisotropic in more itinerant compounds [14] and for 90° 
bond angles [10]. On the honeycomb lattice of Na2lr03, 
for example, Ising terms couple a different spin compo- 
nent S x , S y , or S z along the three equivalent lattice di- 
rections, as in the famous 'Kitaev' model [2, 15] with its 
anyonic excitations. 

In this Letter, we show that such a 'Kitaev'-type 
anisotropy on the triangular lattice, a geometry very 
close the the honeycomb case Na2lr03, induces an in- 
commensurate (IC) and non-coplanar spin pattern that 
turns out to be a lattice of Z2 vortices. The vortices are 
topologically stable point defects of SO(3) order param- 
eters [17], as have long been known to be thermally in- 
duced in frustrated Heisenberg antiferromagnets [17-21] 
as well as in the superfluid A-phase of 3 He [22, 23], and 
more recently in spinor Bose-Einstein condensates [24]. 
In contrast to Z vortices driving Berezinsky-Kosterlitz- 
Thoulcss transitions, clockwise and counter-clockwise ro- 
tation of a Z2 vortex are not topologically distinct, and 
two Z 2 vortices can thus annihilate each other. 




FIG. 1. (color online) (a) Triangular lattice with three types 
of bonds, "xx", "yy" and "zz", see Eq. (1). Sites 1-4 specify 
the four sublattices for the spin rotations S o S [13]. S is de- 
fined as Si = Si,_S a =_(-§$, -S y 2 ,S z 2 ), S 3 = (-Sf , Sf , -SI), 
and S4 = (S4, —5*1, — S4). See (b) for an illustration of S cor- 
responding to S || (1, 1, 1). (c) Spins illustrate two states 
in the "frustration-protected" compass manifold around the 
AFM Kitaev point ip — n/2: Horizontal AFM chains arise 
through coupling oc Jk- For any coupling \Jh\ <C Jk along 
xx and zz bonds, contributions from up and down spins cancel 
and either orientation of the chain has the same energy [16]. 
Darker shading designates the dual lattice of every third up- 
ward pointing triangle, the vorticity in Fig. 3(e) was obtained 
for Wilson loops like the rhombus [17]. 

Here, however, the vortices are not induced by tem- 
perature T but persist for T — ¥ 0. Their density is de- 
termined by the anisotropic couplings and they form a 
regular lattice maximizing their distance. We thus con- 
clude that creating a vortex lowers the energy and that 
the vortices repel each other. The spin-structure factor 
is strongly peaked at three IC momenta, but has addi- 
tional signals. The state may be related to an enigmatic 
phase of a strongly correlated quantum-spin-Hall (QSH) 
model [25]. 

We study the Hamiltonian 

H=J H J2^i-S j + J K Yl S ? S J (!) 

7)( < iJ)IIt 

where Jr is an isotropic Heisenberg coupling between 
classical spins Si on nearest neighbor (NN) sites on 
the triangular lattice. The terms oc Jk are Ising-like and 
couple a different spin component S 7 (S x , S v , and S 2 ) on 
bonds along the three lattice directions 7, see Fig. 1(a). 
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FIG. 2. (color online) Classical phase diagram of Eq. (1). The 
inset hexagons denote the first Brillouin zone, with colored 
markers (resp. line segments at ip = ±7r/2) indicating the 
minima of the Fourier-transformed interactions i7 7 (Q) [31]. 
Arrows indicate how they move when going around the circle 
in a counter-clockwise direction. The lighter shading around 
ip = it/ 2, within the incommensurate regime, indicates a 
phase where Monte Carlo simulations show nematic correla- 
tions over wide temperature ranges. The right and left sides 
can be mapped onto each other by the four-sublattice trans- 
formation [13], mapping FMf>FM and / <H> /. The Kitaev 
points ip = ±7r/2 map to themselves. Stars designate the hid- 
den SO(3) points AFM and FM at tan^ = -2. The AFM 

and AFM points are surrounded by incommensurate phases 
/ and I. 



With Jk = —2Jh > 0, Eq. (1) was proposed for cobal- 
tates [13] and also arises in a strong-coupling limit [26, 27] 
of a QSH model [25, 28], see the discussion toward the 
end of this Letter. For iridium ions in octahedral coor- 
dination and the bond geometry of the triangular lat- 
tice, ferromagnetic (FM) Jh < and antiferromagnetic 
(AFM) Jjc > have been argued to arise from virtual 
hopping to e g orbitals [29] , while other processes lead to 
Jh > and Jk < [10]. As cither sign can thus be 
argued to be relevant for either coupling, we investigate 
all relative strengths by parametrizing Jh = cos ip and 
Jk = sin ip, with ip g [0, 2tt). In real space, Eq. (1) has 
translational invariance plus a two-fold axis perpendic- 
ular to the lattice. In spin space, SO (3) symmetry is 
explicitly broken by the Kitaev terms and only mirror 
planes remain. Combining real and spin space gives the 
D 3d (3m) point group, the high-T space group of many 
layered triangular compounds of ABO2 type. 

First information about the classical phase diagram 
can be obtained by minimizing the Fourier transforms 
J^(Q) of the interaction [30, 31]. As can be seen in 
Fig. 2, the minima for the three spin components 7 can 
be located at different IC momenta. Before going to the 
IC regime, which hosts the Z2-vortex phase, we discuss 
special points and commensurate phases. 



(Hidden) FM and AFM SO (3) symmetric points.— 
The pure Heisenberg points for ip = 0, 7T (Jk = 0) have 
SO(3) symmetry. The FM state at ip — tt is moreover 
an eigenstate of the full quantum-mechanical model; the 
120° pattern [32] at the AFM Heisenberg point ip = is 
not, but is expected to be stable against quantum fluc- 
tuations at T = [33]. 

For tan?/' = —2, the system has a "hidden" SO(3) 
symmetry that becomes apparent after a four-sublattice 
spin rotation [13], see Fig. 1. The transformation pre- 
serves the form of the Hamiltonian (1), but changes Jh —> 
J' h = -Jh and J k -> J' k = 2J h + J k . The AFM point 
ip ~ 0.64767T (where J' h > and J' k = 0) was discussed in 
Ref. [13] and more recently in Refs. [26, 27]: The three 
sublattices of the 120° state together with the four sub- 
lattices of the transformation S f-> S give a twelve-site 
unit cell with generally non-coplanar spins [34] . 

The FM counter point at ip ~ — 0.352-7T has unit cells 
containing up to four spins, with AFM spin patterns de- 
pending on the direction chosen by the rotated spin S. 
The FM and FM phases are quite stable against devi- 
ations from their respective (hidden) SO(3)-symmetric 
point, see Fig. 2. As on the honeycomb lattice [15], quan- 
tum fluctuations favor a collinear AFM phase [27]. 

Kitaev points and nematic ordering at finite T. — For 
ip = tt/2 (3tt/2) and Jh — 0, the lowest energy E = 
— \Jk\ is reached by forming AFM (FM) Ising chains 
along one of the three lattice directions, with spins point- 
ing along x, y, or z [see 1(c)]. Flipping spins of a chain 
does not change the energy, leading to a degeneracy of 
3 x 2 L states growing with the linear system size L. The 
sub-extensive degeneracy is related to a symmetry that 
is intermediate between local and global [35]. An acci- 
dental degeneracy actually allows further ground states 
with spins canted away from the axes [36], but Montc- 
Carlo (MC) simulations show this accidental degeneracy 
to be lifted by thermal fluctuations [31]. Such partial or- 
der into decoupled chains and a resulting nematic phase 
can also emerge in the low-energy limit of more com- 
plex models [37] and are in agreement with the analogous 
square-lattice 'compass' model, as is the entropic lifting 
of an accidental degeneracy [38, 39]. 

On the square lattice, two-dimensional long-range Neel 
order easily replaces the decoupled chains when the 
Hamiltonian deviates from the exact compass form [40, 
41]. This remains true for our FM Kitaev point Jk = — 1 
(ip = 37r/2): A small Jh ^ couples the FM chains with 
an energy oc JhL and stabilizes the FM (FM ) state for 
Jh < (Jh > 0), leading to the sharp first-order transi- 
tion between the FM and FM phases seen in Fig. 2. 

Around the AFM Kitaev point J K sa 1 (ip rj tt/2), 
on the other hand, each spin's neighbors in an adja- 
cent Ising- AFM chain compensate each other [16], see 
Fig. 1(c). MC simulations for ip = atan 4, ip = atan 10, 
and ip = atan 20 support this, as they show partial order 
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FIG. 3. (color online) (a) Spin pattern obtained by MC (down 
to ft = 96) and subsequent optimization for tanip — —0.3 
and 48 x 48 sites, (b) Spin pattern on every third site and 
(c) spin-structure factors 5 7 (k) of (a), (d) shows the vector 
chirality rc(r), Eq. (3), where vortices can be seen; (e) the 
corresponding vorticity of n(r) indicates strongly localized 
Z2 vortices [31]. (f) gives the local energy e Ti — Jh S Ti ■ 

S Ti+ 6 + Jl< S\ 



Sf.ig, where 8 denotes NN bonds. In 



(a), (b), and (d) shading refers to the out-of plane component 
from —1 (dark) to 1 (light yellow); for visibility, only part of 
the lattice is shown here, see [31] for full plots. 



into decoupled chains over wide temperature ranges. At 
very low T, however, spins cant slightly away from the 
preferred direction and an IC non-coplanar ground state 
is selected. 

Z2-vortex phase. — The most intriguing phase is found 
around the AFM and AFM points. As they can be 
mapped onto each other [13], we focus here on the for- 
mer, and in particular on ip < (Jh > 0, Jk < 0), which 
corresponds to the scenario of Ref. 10. For the 120° state 
at ip = 0, the spin-structure factors 



5T(k) Y, , '' k "' ' S 



7|2 



(2) 



are all peaked at the two K points, the corners of the 
first Brillouin zone. This state is also characterized by 
|/c(r)| = 1, where 
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K(r) = —= (S r X S r+ai +S r+ai X S r+a2 +S r+a2 X S r ) , (3) 

3v 3 



is the vector chirality on upward pointing triangles 
[shaded in Fig. 1(c)], with primitive lattice vectors ai^. 

When FM Kitaev-coupling Jk is switched on, spin pat- 
terns remain locally close to the 120° pattern, but are 
distorted at larger distances, see the spin configuration 
in Fig. 3(a) for tani/j = —0.3, which was obtained by MC 
down to ft = 96 and subsequent spin optimization. The 
incommensurate modulation becomes clearer when look- 
ing at only one of the three sublattices of the 120° state, 
see Fig. 3(b). In momentum space, the dominant peaks 
of iS 7 (k) each move slightly away from K, see Fig. 3(c). 
Additional harmonics are small, but do not decrease for 
T — >• 0. These additional features, which have been sug- 
gested to assist the formation of more complex states over 
spirals [42], allow the spin length to be 1 at every site 
and can be taken as a hint in favor of localized defects. 
Indeed, defects are clearly revealed in the site-dependent 
energy shown in Fig. 3(f) as well as in \n(r)\ (not shown), 
which is reduced from 1 to ps 0.85 at the same locations. 

Since we are close to the 120° state, whose order pa- 
rameter has an SO (3) space with the first homotopy 
group 7Ti(SO(3)) = Z2, localized topological defects are 
expected to be Z2 vortices [17]. Following Ref. 17, we 
thus calculate the rotation of n(r) when going around 
loops like the thick blue rhombus in Fig. 1(c) [31]. As 
Fig. 3(e) shows, the vorticity vanishes for most plaque- 
ttes, or adds to one on two neighboring plaquettes, indi- 
cating a rotation by 2-7T [43] . 

The formation of Z2 vortices is further elucidated by 
the form of the continuum limit. Following Ref. [44], a 
classical action for the long wavelength limit of Eq. (1) 
can be derived. We find that the Kitaev terms gen- 
erate Lifshitz invariants, i.e. terms linear in gradients 
of the order parameter [45], in addition to exchange 
anisotropy terms. The gradient terms involve all three 
lattice directions and should dramatically alter the long- 
distance behaviour as compared to the pure AFM model 
(Jj£ = 0). The continuum limit implies that potential 
solutions are either (i) helicoidal one-dimensional modu- 
lations usually found in frustrated systems with Lifshitz 
invariants [45], or (ii) double-twisted states with local- 
ized cores as known from other highly symmetric sys- 
tems, where Lifshitz invariants involve gradients in sev- 
eral spatial directions, as the blue phases of cholesteric 
liquid crystals [46] and chiral magnets [47]. Condensates 
of double-twisted states can create topological defects, 
which fits with the observation that the ground state of 
the lattice model (1) is given by Z 2 vortices condensed 
into a defect-ordered phase, instead of a coplanar one- 
dimensional spiral state [48]. This suggests that the 
'double- twisting mechanism' is operative around the Z2 
vortices, and that the energy gain due to the twisting is 
more important than the loss at the vortex core, where 
the lattice cut-off possibly helps in bounding the latter. 

Around the 120° state at ip = 0, MC simulations con- 
sistently yield the Z2 phase down to the lowest T simu- 
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lated, both for ip < and ip > 0. However, for some 
< ^ 7r /4, we found vortex- free states with very 
slightly lower energy on smaller system sizes L < 48, e.g., 
-1.588 vs. -1.585 for ip = 0.15vr = atan (0.5). Around 
the AFM Kitaev point ip = tt/2, the ground state remains 
IC and non-coplanar, but one spin component 7 domi- 
nates. At low to intermediate T, a 'frustration-protected' 
manifold of (almost) decoupled AFM chains induces a 
nematic phase, see above and Fig. 1(c). Implementation 
and observation of ' compass- model'-type physics, e.g. for 
the realization of protected qbits [8, 9, 41, 49, 50], might 
thus be more promising on such a frustrated than on a 
square lattice. It has not yet been clarified whether the 
transition between the two different IC phases is con- 
tinuous or not, its approximate location is indicated by 
shading in Fig. 2. 

Correlated QSH model. — As mentioned in the be- 
ginning, a situation similar to our model arises in the 
strongly correlated limit of a model for honeycomb iri- 
dates [25, 28]: Next-nearest neighbor sites, which form 
two triangular sublattices, are coupled by a combination 
of Ising and Heisenberg terms oc J2 that bring each sub- 
lattice close to the AFM point [26, 27]. AFM Heisenberg 
NN exchange J\ couples the two sublattices and medi- 
ates an effective FM interaction within each, because the 
Neel AFM phase corresponds to FM order within a sub- 
lattice. When going from dominant Jj to the AFM phase 
at large J2, the functional renormalization group has re- 
cently revealed IC magnetic order [26] for the quantum 
model. We have verified that the situation within the 
sublattices is indeed very close to our triangular lattice, 
in particular, the IC phase of the classical honeycomb 
model is likewise non-coplanar and Z2-vortices can be 
identified. 

Conclusions. — We have investigated an anisotropic 
spin-orbit model given by Kitaev-like couplings on a tri- 
angular lattice, which also describes some features of cor- 
related spin-orbit-coupled electrons on a honeycomb lat- 
tice. In addition to nematic and commensurate phases, 
the classical limit of the model has IC and noncoplanar 
phases characterized by localized defects, topologically 
stable Z2 vortices. In contrast to well-known entropy- 
driven Z2 vortices [17, 23], the defects are here stabilized 
by the 'Kitaev'-type anisotropic spin-orbit couplings and 
persist for T — > 0. We have thus found a Z2 analog to 
the Z-vortex phases of type-II superconductors [52] and 
liquid crystals (the twist-grain-boundary state [53]), and 
to the Skyrmion lattices in chiral magnets [54]. 

Such localizcd-defect states could be realized in tri- 
angular quasi-two-dimensional lattices of 5c? ions (e.g. 
iridium) with strong spin-orbit interaction and electronic 
correlations. It is noteworthy that the zigzag pat- 
tern of Na 2 Ir0 3 [51] is the AFM phase of the Kitaev- 
Heisenberg model on the honeycomb lattice [29] , the cou- 
pling regime of potential triangular-lattice compounds 



can consequently be expected to be close to the AFM 
point as well, i.e., within the regime of the Z2-vortex 
lattice. Optical lattices modelling spin-orbit coupled 
systems might similarly provide a route [55]. While 
we have concentrated on a spin(-orbit) model, itiner- 
ant carriers raise further possibilities, because a non- 
coplanar spin-background can induce topologically non- 
trivial bands [56, 57]. It is thus intriguing to consider 
relations between the Z 2 -vortex phase on the honeycomb 
model to an enigmatic QSH* state of the corresponding 
Hubbard model [25]. 
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Supplemental material 




FIG. 4. (color online) Full-lattice versions for Figs. 3(a), (b), 
and (d) of the main text. 

In this supplemental material, we provide some techni- 
cal notes and additional data supporting the conclusions 
presented in the main text. As our most prominent point 
is the discovery of a Z2-vortex lattice, the first Section 
gives details on how the Z2 vortices were identified from 
the Markov-chain Monte-Carlo (MCMC) data. The sec- 
ond part discusses issues related to the other phases and 
the phase diagram, Fig. 2 of the main text. 

Identification of Z2 vortices 

Figure 4 shows the whole 48 x 48 site lattice for the 
figures shown in Figs. 3(a), (b), and (d) of the main 
text, where only part of the lattice is displayed for vis- 
ibility. The vorticity of n(r) was calculated by going 




FIG. 5. (Color online) Upward pointing triangles (shaded) of 
the triangular lattice, for which «(r) is calculated. Dark shad- 
ing denotes every third upward-pointing triangle, the vorticity 
of #t(r) is evaluated on loops connecting these darker trian- 
gles. The smallest possible loops are right- and left-pointing 
triangles, information from both has to be combined to cover 
the original lattice. When using rhombus loops connecting 
four dark triangles, one rhombus-lattice covers the whole lat- 
tice. These latter loops were used to obtain the figures shown 
here and in the main text. Arrows illustrate the 120° pattern, 
the vector chirality fc(r) on shaded triangles points out of the 
plane. 

around loops as shown in Fig. 1(c) of the main text, i.e., 
from k(t) on every third of the upward-pointing trian- 
gles. This dual lattice in turn forms again a triangular 
lattice, the geometry of the original and dual lattices can 
be seen in Fig. 5. 

We closely follow Ref. 1, with one exception: our vor- 
ticity is rescaled to be for no rotation and 1 for a 
rotation by 2tt, as opposed to 1 and —1, resp. In or- 
der to determine the rotation between two plaquettes of 
the dual lattice connected by an edge of the loop, we 
need two vectors on each plaquette, as a single vector 
«(r) would not determine the rotation axis. We tried 
three approaches: (i) using one spin on each plaquette 
as the second vector, [1] (ii) using /c(r) of the adjacent 
downward pointing triangle, and (iii) not using a second 
vector and instead choosing the rotation axis as orthog- 
onal to the plane spanned by the two K(r)-vectors. All 
three approaches lead to equivalent results near the 120° 
order, and all three suffer from difficulties once the vor- 
tices become very dense and a local 120° pattern can 
no longer be robustly assigned. Figures presented were 
obtained with approach (iii). 

We additionally calculated the vorticity for smaller 
loops connecting only three plaquettes that form right- 
and left-pointing triangles, see Fig. 5. The results are 
equivalent. We also verified that the results remain un- 
changed, when the dual lattice is constructed from the 
downward pointing triangles of the original lattice. 

Fourier transform of the Kitaev-Heisenberg 
interactions 

Since we are dealing with a Bravais lattice and classical 
spins, a promising start is the so-called Luttinger-Tisza 
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method. [2] With Sk = ^ r elkr S r , the classical energy 
per site reads 

e d = 5^7( k ) S k 5 ^ki where (4) 

k.7 

J x (k) = J h [cos(k-l g ) +cos(k-l 2 )] + (J h + J fc )cos(k-l x ), 
J y {k) = J x (R^-k), J z {k) = J x {R_^-k) . (5) 

Here l y = x, 1 2 = (x + v3y)/2, l x = l z ~ l y are primitive 
vectors of the lattice, and R ± 2w stand for C3 rotations 
around the z-axis in momentum space. For a large regime 
of — atan(3/2) < ip<3ir/2, all three i/ 7 (k) are minimized 
at k = (0,0), indicating a FM state. Here the ground 
state (GS) energy does not change by global rotations 
of the spins, and so the GS manifold has a continuous 
accidental degeneracy (except for ip = 7r where the de- 
generacy is related to the global SO (3) symmetry). 

For the range 37r/2 < ip < — atan(l/2) the three minima 
are at the M points of the first Brillouin zone, i.e., the 
midpoint of the sides of the hexagon, with each J7 7 being 
minimum at a different side. This is the signature of the 
commensurate order that can be seen as the FM phase. 
Since the three M points are one-half of a reciprocal vec- 
tor, one can generate a continuous family of ground states 
by forming linear combinations of the three modulations 
without violating the unit length constraint. [3] The re- 
sulting GS manifold has a two-dimensional structure con- 
tinuously connecting the three M points. In the rotated 
frame (defined by the transformation S — > S, see main 
text), this is nothing but the accidental degeneracy re- 
lated to global rotations (again, except for the FM point 
where the degeneracy is symmetry related). 

Directly at the pure Kitaev points ip = tt/2,3tt/2, each 
Sf-yik) is minimal along a line in momentum space, a 
signature of the decoupled chains and the "intermediate 
symmetry" (see main text), as in the square-lattice com- 
pass model. The GS manifold has an additional contin- 
uous accidental degeneracy which will be discussed else- 
where. 

In the remaining range — atan(l/2) < ip < — atan (3/2), 
each J-yik) is peaked at a different incommensurate [4] 
momentum k = Q 7 for each spin component 7, which 
are related by C3 rotations. Consequently, not all three 
spin components can be minimized simultaneously with 
a single Q 7 . If one were able to do so, or even if one 
were able to optimize only two components 7 at the same 



k = Q, the problem could be solved by a coplanar spiral 
S r /5 = cos(Q-r)ei + sin(Q-r)e2, where the two vectors 
ei and e2 define the spin plane. As there is however no 
momentum that minimizes simultaneously the energy for 
more than one spin component, there is no way to form a 
ground state through a coplanar incommensurate spiral. 
Compass-like phases at the Kitaev point 

Figure 6 shows that - as in the square-lattice compass 
model - low-T states have spins polarized along one of the 
three axes, even though the ground-state manifold has an 
additional accidental degeneracy and contains also non- 
collincar and even non-coplanar states. 
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FIG. 6. (color online) Spin-structure factors <S 7 (k) for ip = 
7r/2 and inverse temperature (a) /3 = 1/fcsT = 8, where 
signals in all three components are visible, (b) pi — 10, (c) 
P = 100, and (d) /3 = 1000, where decoupled chains along one 
direction dominate. 



